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Abstract 

The dynamics of Bose-Einstein condensates trapped in a deep optical lattice 
is governed by a discrete nonlinear equation (DNL). Its degree of nonlinearity 
and the intersite hopping rates are retrieved from a nonlinear tight-binding 
approximation taking into account the effective dimensionality of each conden- 
sate. We derive analytically the Bloch and the Bogoliubov excitation spectra, 
and the velocity of sound waves emitted by a traveling condensate. Within 
a Lagrangian formalism, we obtain Newtonian-like equations of motion of 
localized wavepackets. We calculate the ground-state atomic distribution in 
the presence of an harmonic confining potential, and the frequencies of small 
amplitude dipole and quadrupole oscillations. We finally quantize the DNL, 
recovering an extended Bose-Hubbard model. 
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I. INTRODUCTION 



Bose-Einstein condensates (BECs) trapped in periodic potentials are a new bedtest to 
investigate different and fundamental issues of quantum mechanics, ranging from quantum 
phase transitions [1] and atom optics [2,3] to the dynamics of Bloch and Josephson oscilla- 
tions [4-6]. 

The key feature of these systems is the competition/interplay between the discrete trans- 
lational invariance (introduced by the periodic potential) and the nonlinearity (due to the 
interatomic interaction). For instance, it has been pointed out that the excitation spectrum 
exhibits a band structure which has analogies with the electron Bloch bands in metals [7-9] . 
However, the coexistence of Bloch bands and nonlinearity allows for solitonic structures [10] 
and dynamical instabilities [11-13] which do not have an analog neither in metals, nor in 
Galilean invariant nonlinear systems. 

The BEG in a periodic potential is described in mean field (or classical) approximation 
by the Gross-Pitaevskii equation (GPE) Eq.(4). When the interwell barriers of the periodic 
structure are high enough (see below) , the relevant observables of the system are the number 
of particles Nj{t) and the relative phases ipj+i{t) — (pj{t) of the condensates in the lattice 
(the subscript j denotes the different wells of the array). In [10] it has been shown that 
the amplitudes ipj = y^iV^e*"^^' satisfy a discrete nonlinear Schrodinger equation (DNLS), 
recovered from the GPE (4) in the tight-binding approximation. Such a mapping has allowed 
the investigation of localized and extended excitations [10,14,15] in the framework of the 
nonlinear lattice theory [16], and has clarified the connection between such systems and 
an array of superconducting Josephson junctions. The relative phase dynamics of different 
condensates has been experimentally investigated looking at the interference patterns created 
by the atoms tunneling in the continuum from a vertical lattice [4], or by the overlapping 
condensates once the confining traps are removed [3,6,17]. 

In the "standard" tight-binding approximation, the condensate wavefunction localized 
in the jth well is factorized as a dynamical part ipj{t) = J Nj{t)e^'^^^*^ and a spatial, real 



2 



wavefunction $j(r) centered in the minimum fj of the well. A crucial assumption, which 
will be discussed in detail, is that the shape of does not depend on the number of 

particles Nj{t) in the same well. Under such conditions, the condensate order parameter 
^{f, t) can be written as: 

^{f,t)^Y.^,{t)^,{r). (1) 
i 

The DNLS equation is recovered replacing Eq.(l) in the GPE Eq.(4), and integrating out 
the spatial degrees of freedom. Neglecting higher order terms, we get [10]: 

ih-^ = -K{tpj-1 + ^j+l) + C/2 I r ^3 + ^3^3 (2) 

with K, U2, ej depending on the geometry of the trapping potentials and on the average 
number of atoms in each well [cfr. with Eqs.(8-ll)]. 

II. A GENERALIZED TIGHT-BINDING APPROXIMATION 

In this paper we stress the importance (and exploit the consequences) to generalize 
the "standard" tight-binding approximation Eq.(l). This generahzation is imposed by the 
nonlinearity of the GPE Eq.(4), and largely extend the range of validity of the DNLS Eq.(2) 
in the study of the dynamics of weakly coupled BECs. The central argument is that the 
density profile of each condensate can strongly depend on the number of atoms present at a 
given instant in the same well. This introduce site- and time- dependent parameters in the 
DNLS Eq.(2), modifying, in particular, its effective degree of nonlinearity. Therefore, the 
tight-binding approximation of nonlinear systems has to be generalized as: 

W,t)^T.Mt)^3ir;N,{t)). (3) 

3 

with ^j{f; Nj{t)) depending implicitly on time through Nj{t) = \ipj{t)\'^. We stress here, 
and discuss again later, that the spatial wavef unctions $j (which are considered sufficiently 
localized in each well) can also depend explicitly on time due to the excitation of internal 
modes. In the present approach, however, we consider the adiabatic limit in which the 
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interwell number/phase dynamics is much slower that the typical time associated with the 
excitations of such internal modes (and, of course, the cases where such modes are not 
already present in the initial configuration of the system). In this limit, which can is well 
satisfied in typical experiments, the spatial wavefunctions in Eq.(3) will adiabatically follow 
the tunneling dynamics and can be approximate with the rea/ wavefunction ^j{'r; Nj{t)). 

III. THE DISCRETE NONLINEAR EQUATION 

The BEC dynamics at T = satisfies the Gross-Pitaevskn equation (GPE) [18]: 

ih^ = -^^'"^ + [Ve.t + ^0 I * ri* (4) 

where Vext is the external potential and go ~ ° , with m the atomic mass and a the 
s-wave scattering length: a > (a < 0) corresponds to an effective interatomic repulsion 
(attraction). For the sake of clarity, we will focus on the case a > (as for ^'^Rb atoms), 
and with Vext will be given by the optical periodic potential Vp superimposed to a harmonic 
magnetic field Vm- The periodic potential is 

Vp^Vosin^ikx) (5) 

where k = 27r/A and A is the wavelength of the lasers (the lattice spacing is A/2). The 
energy barrier between adjacent sites, Vq — sEr, is expressed in units of the recoil energy 
Ep = From (5) we see that the minima of the laser potential are located at the points 
Xj = j |. Around these points, Vp a; ^u;l{x — Xj)"^, where 



2Vok^ ... 
V m 

The magnetic potential is Vm = ^[^Ix^ + ^yU"^ + ^1^"^]) with Cjx » ujx- It is convenient 
to write the external potential as V^xt — Vl + Vd, with the confining lattice potential 
Vl = Vosin'^ikx) + f[a;^?/^ + a;^^^], and the "driving" field Vd = f o;^ x^. Vd has a 
simple physical meaning: F = — is the effective force acting on the center of mass of a 
condensate wave packet moving in the periodic potential, see Section VII. 



Here we consider a one-dimensional optical lattice superimposed to an harmonic driving 
field, but the following considerations can be easily generalized to arbitrary Vd and, in 
particular, extended to the case of two- [19] and three- [1] dimensional arrays created by 
several counterpropagating laser beams. 

Replacing the nonlinear tight-binding approximation (3) in the GPE (4) and integrating 
out the spatial degrees of freedom we find the following DNL: 



(/lb ■ 

-[K+x (I V. r + 1 T)] i^j+i - [K+x (I r + 1 n V'.-i (7) 



(we use the normalization / $^ = 1, while the total number of atoms is Nt — J2j iV'jP = 
J2j Nj)- In Eq.(7), the "local" chemical potential is the sum of three contributions 

^toc (Jepends on the atom number A^^ through the condensed wavefunction $j(r; Nj{t)). The 
tunneling rates Kjj±i between the adjacent sites j and j ± 1 also depend, in principle, on 
the respective populations: Kjj±i{Nj; Nj±i) = — J df (^V$j • V%±i + $jVea;t^i±i)- In 
this case, however, we can expand the wavefunctions around an average number of atoms 
per site, A^o, and keep only the zero order term ^j{Nj) ~ ^j{No): 



I '^^ ^2;^^^^ ■ ^^^^^ + ^.Kxi^.-ii]. (9) 
We have checked numerically that higher order terms are negligible: for instance, with the 

dK 
dNo' 



K 



experimental setup of [6], Vq = IOEr and A^o = 1000, we have Ki = if^N ~ IQ-^X. 



Similarly, the coefficient x is given by: 

x^-gJ rfr $^$,±1. (10) 



The on-site energies arising from any external potential superimposed to the optical lattice 
are 

e, = j dfVD (11) 



Cj oc {sj oc j) when the driving field is harmonic (hnear). We notice that in the hmit 
ujx » u)x considered here, ej does not depend on the on-site atomic populations. 

In the derivation of Eq.(7) we have exploited the (quasi-)orthogonality of the con- 
densate wavefunctions / df $j ^j±i — 0. Moreover, we have verified numerically that 
spatial integrals involving condensates distant more than one site, as well as terms pro- 
portional to / df $1 can be neglected. E.g., with Vq = ^OEr and Nq = 1000, 
goNoJ $2 ^^.^JK ~ 10-^ while x^o/K ~ IQ-^. For Vo = 20Er and A^o = 10000, 
xNq/K ~ 10-^ In a double well with e.g. o;^ = 27r • lOOHz and A^o = 10000, x^o ~ K, 
while Ki and go Nq J df $| ^j±i can still be ignored. For these reasons, we do not neglect 
the X terms in Eq.(7). A detailed account of the related numerical study will be presented 
elsewhere. 

The atom number dependence in Eq.(8) introduces an effective time- dependent, real, 
local chemical potential ijL''j"^[Nj{t)]. This refiects an important approximation contained 
in the DNL: terms proportional to d^j/dt have been neglected. In other words, we have 
neglected the phases associated with the spatial dynamics of the $j(r; Nj{t)) in Eq.(3). This 
adiabatic approximation is well satisfied when the tunneling time (~ Nj/Nj) is much longer 
than time scales associated with the change in shape of the wavefunctions (~ uj~^, oj~^). 
In this limit, well satisfied in realistic experiments, the spatial profile of the wavefunctions 
adapts adiabatically to the instantaneous number of atoms present in the respective wells 
[20]. 

The dependence of the local chemical potential on the number of atoms depends on the 
effective dimensionality of the condensates trapped in each well of the lattice. This can 
be determined comparing the interaction chemical potential //^"* = iV'iPfi'o / and the 

three frequencies, Ux, ^y, obtained expanding the lattice potential around the minima 
of each well Vl — y [a)^(a; — XjY -\- Uyu"^ -\- ojIz^]. A sufficiently accurate calculation of ji^* 
as a function of Nj can be obtained approximating the condensate order parameters with 
gaussians or Thomas- Fermi functions [21]. Here we first consider some limit cases which are 
particularly instructive. 
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When ojj^^ijjy^LiJz » the spatial widths of each trapped condensate do not de- 

pend (in first approximation) on the number of particles Nj in the same well, and the 
condensates wavefunctions are well approximated by gaussians. We consider this as a OD 
(zero-dimensional) case {nD, with n — 0,1, 2, 3, should not be confused with the spatial 
dimensionahty of the lattice), and the ansatz Eq.(3) reduces to the ordinary TBA Eq.(l). 
The ID case arises when two frequencies are greater than the interaction chemical potential. 
For instance, \iLijx,0Jz » A*}"* » ^y, the system realizes an array of weakly coupled cigar- 
shaped condensates oriented along the y-sods: the wavefunction $j will be factorized as a 
product of two gaussians (in the x and z directions) and a Thomas-Fermi in the y variable. 
In the 2D case only one frequency is smaller than the local interaction chemical potential. 
If (Da; >> >> u}y,uiz-i we have an array of pancake- like condensates, with <I>j factorized 
as a gaussian (along x) and a Thomas- Fermi in the y and z variables [see Eqs. (31-33)]. The 
3D case is given by the condition //^"* >> Lijx,ojy,uJz and the wavefunction in the jth well 
$j is simply given by a three-dimensional Thomas-Fermi function. To summarize: 
2>D ij,"^* {r^ Nj^^) » u)a,0Jb,uJc [spherical] 



2D Ua » Atf* (~ N]''^) » ujb,ujc \pancake] 



ID uja,ujb » /xf * (~ Nf^) » Uc [cigar] 
OD uja,uJb:<-^c >> IJ'^j^^ Nj) [spherical] 
with a,h,c arbitrarly corresponding to the x,y,z spatial directions, and among square brack- 
ets are specified the geometric shapes of the condensates in each well. 

The crucial point is that the effective dimensionality of the condensates gives a different 
scaling of the local interaction chemical potential Eq.(8) with the number of atoms: 

= I r (12) 

a = ^^, D = 0,1,2,3. (13) 

where Ua is a constant which does not depend on the number of atoms nor on the site 
index. In the following we will often consider, for the sake of clarity, the limit cases when 



the local chemical potential is given by Eq.(12) (generalization to more complicate functional 
dependences of fi^J"^ from A^^ are straightforward). The DNLS Eq.(2) is recovered from the 
DNL Eq.(7) in the case D = (i.e. a = 2) and neglecting terms proportional to x- 

The derivation of the Hamiltonian of the system requires some care. The dynamical 
variables t/jpiht/jj are canonically conjugate (■^j = d(ihi*) ) "^i^^ respect to the effective 
Hamiltonian 

^e// = Y.\ ej^jipj - K{il:*'il:j+i + c.c.) 

3 ^ 

-x[ I r + r^-i) + c.c] + ^t/^i^^-r+^l (14) 

(with the nonlinear term f^^f^alV^jl"^^ obtained from tpj J dip* [j}-"^) . 

The effective Hamiltonian Ti-eff is an exact integral of motion, but differs from the 
"adiabatic" Hamiltonian retrieved simply replacing Eq.(3) in the Gross-Pitaevskii energy 
functional: 

T-iad = - + c.c.) + 

j ^ 

-X[ I P ^,(^*+i + + C.C.] + ^t/|^,f |, (15) 

with U = goJdr^'j. The ?Ye// and Ti^d are identical only in the OD case. In general, Had 
is not exactly, but only "adiabatically" , conserved during the dynamics. 



IV. EXCITATION SPECTRA 

We now derive the Bloch excitation spectra and the Bogoliubov dispersion relation 
of the system (with ej = 0), calculate the sound velocity and investigate the dynamical 
stabihty of condensate travehng waves. Eigenfunctions of the DNL are the plane waves 
i^n — ■00 e^'^~^^^^\ with chemical potential and energy per site: 

= /x'"^ -2(K + AxNo) cosp 
E = E^''^ -2{K + 2x No)No cosp (16) 
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where Nq =| V'o \^ A*'"' = l^f%i=fo = t^a I V'o |°, and 0"'' = V'o /#>J°1v.=V'o = 2C/« | 
V'o l""*"^ /(q^ + 2) [see Eq.(12)]. Prom Eq.(16) we can recover the group velocity of Bloch 
waves with quasimomentum p: Vg = = 2 [K + 2 x A^o) sin p. 

We remark that the Bloch energy E and the chemical potential ji have the same cosp 
dependence on the quasimomentum p, but with different coefficients. This introduces dif- 
ferent effective masses for the system (see also [23]), which will enter in peculiar ways in 
the equations discussed in this paper, as will be reported elsewhere [22]. Here we will write 
down our results only in terms of the DNL parameters. 

In order to derive the Bogoliubov dispersion relation of the system, we perturb the large 
amplitude wave as -^n = [V'o + u{t)e^'^'^ + v*{t)e~'^'^^] e*(*'"~'^*/^) . Retaining only first order 
terms proportional to it/V'o and v/V'o; we get 

with a = 2{K + 4xA^o) sin p sing, b = 2K cosp — 2{K + 4xA^o) cospcosg + iVo^^ and 
c = — 4x'0o cosp (1 + cos q) + i^l^j^- Up to the order x^N^/K'^, we get the eigenvalues: 

cu — 2{K + 4xA^o) sinpsin g± 

2^AK{K + 8xiVo) cos^psin^ | + 2(X + 2xNo)^No cospsin^ |. (18) 



V. SOUND-WAVES AND INSTABILITIES 

The small q (large wavelength) limit of the Bogoliubov dispersion relation Eq.(18) is 
linear. Therefore, the system support (low amplitude) sound waves (propagating on top of 
the large amphtude travehng wave ■00 e*^^""'**/'^)) having velocity 



Vs = ^|,=o = 2{K + 4xA^o) sinp ±^2{K + 2xiVo)|^iVo cosp (19) 

with // given by Eq.(16). The -|- (— ) sign corresponds to a sound wave propagating in 
the same (opposite) direction of the large amplitude traveling wave. Notice that, contrary 
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to the case of a Galilean invariant system (Vq — 0), the sound velocity depends on the 
quasimomentum p. Moreover, Vg depends on the effective dimensionality of the condensates, 
since (from Eq.(12,13)) ^Nq ^ a N^^^. 

In the limit Vq = 0, the system is energetically unstable if a; < 0, namely when the group 
velocity is larger than the sound velocity (Landau criteria for superfluidity). This instability 
is present also when the the system has a discrete translational invariance (Vq > 0): from the 
Bogoliubov excitation spectrum Eq.(18) and the condition a; < 0, we have that the system 
is not superfluid when 

[2{K + AxNo) sinp]2 > 2{K + 2xN^)^N^ cosp (20) 

There is a further, different (dynamical) instability mechanism, which disappears in the 
translational invariant limit (when a > 0) . This instability is associated with the appearance 
of an imaginary component in the Bogoliubov frequencies: from Eq.(18), this component 
appears if cosp < 0. This reflects on an exponential increase of the amplitude of the 
perturbation modes, with the consequent strong dephasing and energy dissipation of the 
condensate traveling wave. The unstable modes g, for a given quasimomentum p, are given 
by: 

2 V + -ir) I I 2 < dN,'''- 

For a = 2 and x = we recover the standard DNLS results [24,12]. The onset of ener- 
getic and dynamical instabilities with an arbitrary Vq has been studied numerically in [11]. 
Experimental evidences are reported on [25] . A different manifestation of the instability is 
associated with the seff-trapping of a condensate wavepacket at rest in an optical lattice 
[10]. First experimental results are reported in [26]. 

VI. GROUND-STATE ATOMIC DISTRIBUTION 

We now consider a magnetic harmonic potential superimposed to the optical lattice 
Ej — flj^, with fl — ma;^A^/8. For a large nonlinearity, the ground-state atomic distribution 
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can be calculated from the DNL (7) in Thomas- Fermi approximation, i.e. neglecting the 
kinetic terms proportional to K and x with respect to the nonlinear term: 

where the inversion point is — v/D,. Replacing sums with integrals, we get 



u — 



where Cq = 2^/"+^[r(2/Q; + l)Y/r{A/a + 2) is a numerical constant (F is the Gamma func- 
tion). For a = 2 (OD), C2 = 4/3, while for a = 1 {2D), Ci = 16/15. 

VII. NEWTONIAN DYNAMICS AND SMALL AMPLITUDE OSCILLATION 

FREQUENCIES 

We now study the wave-packet dynamics of a BEG in an optical lattice. We resort to 
a variational approach, previously considered in [10]. Here we use a general variational 
wavefunction 

tPj = y^/(^l^)e'f0-«)+40-O^ (24) 

where ^{t) and a{t) are, respectively, the center and the width of the wavepacket, p{t) and 
5{t) their associated momenta and ]C{cr) a normalization factor (such that J2j Nj — Nt)- f 
is a generic function, even in the variable X — {j — ^)/a. E.g., we can choose f{X) = e~^^ or 
f{X) = (1— X^)^/" (with — 1 < X < 1) to describe, respectively, the dynamics of a gaussian 
or a Thomas-Fermi wave packet. With the Lagrangian C — iU J2j '^j'>Pj ~'^eff we can recover 
the equations of motions for the variational parameters qi{t) = ^{t),a{t),p{t), 5{t), given by 
= 1^. With the variational wavefunction Eq.(24), the Lagrangian becomes: 

— = np^-fiad—- Vd{^, a) - U^-^ + ^^^(^i ^) ^osp + -^^^x(^; ^) cosp (25) 

where Voii^a) = dX p{X)e{aX + 0, Ua = 2[/„J^L/[(a + 2)Jf/'+'], Ij{a;6) = 
J dXf{X + l/2a)f{X - l/2c7)e-^^ and S) ^ J dXf{X + l/2a)f{X - l/2a)[f{X + 

11 



1/2(7) + p{X - 1/2(7)] e^'^^^. Furthermore Ti = JdXp{X), I2 = JdXX^p{X) and 
^NL — I dXf°'~^'^{X) are real numbers which depend on the particular choice of /. Prom 
the Lagrangian equation of motion we get the group velocity ^ and the effective force acting 
on the center of mass of the wavepacket: 

,2K^ . ^, 2xNt, 



[—Ij{a; 5) + -^^^x(^; ^)] sinp 



nt = -^ (26) 

The frequency of small amplitude oscillations of the wavepacket driven by an harmonic field 
ej = (which gives Vd{^, a) = + a^^)), is: 

u;,,,-j^[2K+^^) (27) 

where A^o — Nt/^g [27]. Eq.(27) has been calculated in the hmit of a large width a » 1, 
where Xj{a;Q) ~ Xi and X^(cr; 0) ~ 2X3, with — j dXp{X). The same results follow 
from the exact equation of motion for ^ = J2jjNj and p = <fj+i — <fj, with the latter 
assumed equal for each j along the array (and using the fact that J2j ^NjNj^i ~ Nt and 
J2j NjyjNjNj+i ~ Nf). For X = Eq.(27) coincides with the result in [6]. 

To calculate the quadrupole oscillation frequency we need the equation of motion for the 
width (7 and the conjugate momentum S (still with Vd — ^C"^)- 
J2 2KdIj 2xNTdI^ 

M^-^ - -2Q?a + + cos . + ^ - cos . (28) 

The equilibrium position is given hy 5 — 0, & — 0, ^ — and p — 0. Linearizing around the 
equilibrium for the Thomas- Fermi ground-state (22), and after a lenghty calculation, we get 
the frequency of the quadrupole oscillations: 

na{a + A)XNLf^r^ , SXN0X4 



where Z4 = / dXX'^ p{X) [28]. Eq.(29) shows that the quadrupole frequency exphcitely 
depends on the effective dimensionality of the condensates in each well [29]. Collecting 
Eqs.(27) and (29) we get 
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,2 

dip 



(30) 



When xA^o << -f^, -2^2*: — 2 ^±|. In particular, i^guadr/^dip = 3 in the zero-dimensional 
case, = 0. The 2D, x = result, i^luadr/^dip = 5/2, is in agreement with the results of 
[30]. 



VIII. NUMERICAL ESTIMATES 

We now consider a specific example to further clarify the calculation of the DNL coeffi- 
cients Eqs.(8-ll). Considering the experimental apparatus of [6], we put: a;^; = 27r x 9Hz, 
ujy — ujz = ojr — ^-K X 90Hz, X — 795 nm, En/h — 3.6 kHz and Vo/Er from 2 to 15. Prom 
Eq.(6), we obtain 0^/'^'^ — Vs ■ 7.2 kHz. 

Since uj^ » uJri we find /x^*" -|- ^^^'^ ~ ^<i'a:/2. With an average value of atoms in 
each well A^o ~ 1000 and with Vq = 5£'r, we obtain an interaction chemical potential 
H^p^ ~ h ■ 2kHz, which corresponds, according to the table in Section III, to the 2D case. 
The system can be seen as a horizontal pile of pankakes, having a smaller diameter at the 
border of the pile, dense at the center and more dilute at the surface. In this limit we have: 

$,(f, N^{t)) ^ - ^) (31) 

where = ((7-\/7r)~^^^e~*^^~^^^^/^'^^ is a gaussian with width a (we impose / dx{(f)QY — 
J dydz{(f)j'p)'^ = 1). A variational calculation shows that there is a very weak dependence of 
a on Nj] we therefore assume it as site-independent [17,31]: a = ^^^1/4 . Replacing Eq.(31) 
in Eq.(4) and integrating out along the x direction, we obtain an equation for (t>TF{y-i z)'- 

[-^V| + V(^) + QoNA^^YWe = /^rV$l (32) 

with ^0 — fl'o/'\/27r(T; R — {y, z) is the vector expressing the position in the y-z radial plan 
and V{R) — '^LolR^. In Thomas-Fermi approximation (i.e. neglecting the kinetic terms in 
Eq.(32)), we find 
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The inversion point is R'j_ = 2/x*"Yma;^. Replacing Eq.(33) in Eqs.(8) we obtain 



The on-site energies (11) are given by ej = where Q = '^muj'l^^Y . We have neglected 
the kinetic terms ef^^^ ~ ^ / dRiV R(t)x^p) , consistent with the the Thomas- Fermi approxi- 
mation (33). Using Eqs.(34) we get the DNL (7) with D^2{a^l) and [see Eq.(12)] 



The population distribution in the ground state, according to Eq.(22), is given by 

'^.-i^Yi^-P)- (3a) 



Hnv 



^UJ V j 

The inversion point is = and the discrete chemical potential (23) is is = 
{IbNrU^V^/Wy/^. Therefore 

where d — X/2, a^o — ^Jfi/muj and uj — {ujluj^yf^. The D — 2 ground state Eqs. (36-37) is 
in agreement with [17], previously calculated with a different approach. 

IX. QUANTUM CASE: AN EXTENDED BOSE-HUBBARD MODEL 

The quantization of the DNL equation requires some care. The quantum equation for 
the bosonic gas in an external potential is: 

i;i^^(f, t) = [T + Ve,t + 90^^^^, (38) 

The Gross-Pitaevskii equation (4) can be retrieved introducing the classical field ^ — {^) 

and with (^t^^) ~ (^t)(^)(^). 

In tight-binding approximation 

14 



*(r,i) = ^^,(i)$,(r-) (39) 

3 

(with the bosonic number operator), we obtain the Bose- Hubbard model (BHM) [32,33] 

H^Y.{- + h.c.) + ^(V.jV'ji-i-) + e.-V'ji-} (40) 

j ^ 

We now discuss the case in which the locahzed wavefunction $j in the jth well adiabatically 
depends on the average number of particles in the same well: the generalization to the 
quantum case of (3) is 

^(f,t) = ^^,(t)%(f;iV,(t)) (41) 

3 

where 

= (^J^,). (42) 

Replacing the ansatz (41) in (38), it is easy to recover the quantum equation of motion for 
bosonic operators ipj. Such equations arc generated, with the standard bosonic commutation 
relations, from the extended Bose- Hubbard Hamiltonian: 

-K{^}^j+, + h.c.) - xi^lMA^U + "Pj-i) + ^-c-lj (43) 

with the parameter K^x^^j^U expressed as in the classical DNL Eq.(7). Notice that the 
extended BHM can be alternatively recovered quantizing the classical adiabatic Hamiltonian 
Had Eq.(15) (and not the effective Hamiltonian Eq.(14)). 

X. CONCLUSIONS 

The Gross-Pitaevskii dynamics of a Bose-Enstein condensate trapped in a deep periodic 
potential can be studied in terms of a discrete, nonlinear equation. This mapping allows 
a clear and intuitive picture of the main dynamical properties of the system, which can 
be calculated analytically We have shown that the slopes of the energy and chemical 
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potential Bloch excitation spectra, with respect to the quasimomentum of the condensate, 
are different. We have calculated the Bogoliubov dispersion relation, and studied the sound- 
wave velocity as a function of i) the effective dimensionality of each condensate, and ii) the 
quasimomentum of the carrier wave. Through a Lagrangian formalism, we have recovered 
Newtonian- like equation of motion of localized wavepackets, and the frequencies of dipole and 
quadrupole small amplitude oscillations. We have finally quantized the discrete nonlinear 
Hamiltonian recovering an extended Boson-Hubbard model. 

Note added in proof: An equation similar to the DNL (7) (with a — 2, and including the 
term proportional to / (if has been derived by Oster, Johansson, and Eriksson [34] 

to describe the amplitude of an electric field in an array of coupled waveguides embedded 
in a material with Kerr nonlinearities. 
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